% The MATLAB source code of TCDE is developed in MATLAB R2019b with 64 bits. 
% Please run this code on MATLAB2019b or later with 64 bits. 
% Please cite this article as: Qiang Yang, Zhuo-Yin Qiao, Peilan Xu, Xin Lin, Xu-Dong Gao, Zi-Jia Wang,Zhen-Yu Lu, Sang-Woon Jeon, Jun Zhang . 
% Triple competitive differential evolution for global numerical optimization.
% Swarm and Evolutionary Computation, https://doi.org/10.1016/j.swevo.2023.101450.
function [Global_score,Global_pos,Convergence_curve] = TCDE(SearchAgents_no,FEsmax,lb,ub,dim,fobj,fnum)
lb=-100;
ub=100;
% Initialize parameter
k=3;
F = 0.3;
CR = 0.8;
H = 6;
Arate=1.3;
NPinit = SearchAgents_no;
NPmin = 4;
counteval = 0;

A = [];
nA = 0;
MF = F * ones(H, 1);
MCR = CR * ones(H, 1);
iM = 1;
% Initialize population and Evaluation
X = lb + (ub - lb) .* rand(dim,SearchAgents_no);
fitness = feval(fobj,X,fnum);
fitness = fitness-fnum*100;
counteval = counteval + SearchAgents_no;
[fitness, fidx] = sort(fitness);
  Convergence_curve(1,1)=fitness(1);
countiter = 2;
X = X(:, fidx);
V = X;
U = X;
Chy = cauchyrnd(0, 0.1, 2*SearchAgents_no);
iChy = 1;
rcount=0;
PR=ones(1,H)/H;
while counteval < FEsmax

    Ri=1:SearchAgents_no;
    Rank=k*(SearchAgents_no-Ri)+1;
    Pr=Rank./sum(Rank);
    p=0.085+0.085*counteval/FEsmax;

    MCR(H)=0.9;
    MF(H)=0.9;
    r=randsample(H,SearchAgents_no,true,PR);
    % Crossover rates
    CR = MCR(r)' + 0.1 * randn(1, SearchAgents_no);
    CR((CR < 0) | (MCR(r)' == -1)) = 0;
    CR(CR > 1) = 1;
    % Scaling factors
    F = zeros(1, SearchAgents_no);
    for i = 1 : SearchAgents_no
        while F(i) <= 0
            F(i) = MF(r(i)) + Chy(iChy);
            iChy = mod(iChy, numel(Chy)) + 1;
        end
    end
    F(F > 1) = 1;

    if counteval<0.6*FEsmax
        F(F>0.7) = 0.7;
    end
    if counteval<0.2*FEsmax
        FW = 0.7*F;
    elseif counteval<0.4*FEsmax
        FW = 0.8*F;
    else
        FW = 1.2*F;
    end
    if counteval<0.25*FEsmax
        CR(CR<0.7) = 0.7;
    elseif counteval<0.5*FEsmax
        CR(CR<0.6) = 0.6;
    end
    % pbest index
    pnum=max(2, round(p * SearchAgents_no));
    pbest = randi(pnum,1,SearchAgents_no);
    % Mutation
    r1 =randsample(SearchAgents_no,SearchAgents_no,true,Pr);
    r2 =randi(SearchAgents_no+nA,1,SearchAgents_no);
    PA=[X A];
    rankindex=randperm(SearchAgents_no);
    for i = 1 :3: SearchAgents_no - 2
                pool(1).value=fitness(rankindex(i));
                pool(1).id=rankindex(i);
                pool(2).value=fitness(rankindex(i+1));
                pool(2).id=rankindex(i+1);
                pool(3).value=fitness(rankindex(i+2));
                pool(3).id=rankindex(i+2);
                [~,index]=sort([pool.value]);
                first=pool(index(1)).id;
                second=pool(index(2)).id;
                third=pool(index(3)).id;
                FFFF=[FW(rankindex(i)) FW(rankindex(i+1)) FW(rankindex(i+2))];
                CCCC=[CR(rankindex(i)) CR(rankindex(i+1)) CR(rankindex(i+2))];
                FFFF=sort(FFFF);
                CCCC=sort(CCCC);
                %3
                V(:, third) = X(:, third) + FFFF(3) .* (X(:, first) - X(:, third)) + FFFF(3) .* (X(:, second) - X(:, third));
                   for j = 1 : dim
                            if V(j, third) < lb
                                V(j, third) = 0.5 * (lb + X(j, third));
                            end
                            if V(j, third) > ub
                                V(j, third) = 0.5 * (ub + X(j, third));
                            end
                    end
                        jrand = randi(dim);
                        for j = 1 : dim
                            if rand <= CCCC(3) || j == jrand
                                U(j, third) = V(j, third);
                            else
                                U(j, third) = X(j, third);
                            end
                        end
                   % 2            
        while second == r1(i)
            if rcount<10
                r1(i) =randsample(SearchAgents_no,1,true,Pr);
                rcount=rcount+1;
            else
                break;
            end
        end
        rcount=0;
  
        while second == r2(i) || r1(i) == r2(i)
            if rcount<10
                r2(i) =randi(SearchAgents_no+nA);
                rcount=rcount+1;
            else
                break;
            end
        end
        rcount=0;
                        V(:, second) = X(:, second) + FFFF(2) .* (X(:, first) - X(:, second)) + FFFF(2) .* (X(:, r1(i)) - PA(:, r2(i)));
                   for j = 1 : dim
                            if V(j, second) < lb
                                V(j, second) = 0.5 * (lb + X(j, second));
                            end
                            if V(j, second) > ub
                                V(j, second) = 0.5 * (ub + X(j, second));
                            end
                    end
                        jrand = randi(dim);
                        for j = 1 : dim
                            if rand <= CCCC(2) || j == jrand
                                U(j, second) = V(j, second);
                            else
                                U(j, second) = X(j, second);
                            end
                        end
                        %1
        while first == r1(i)
            if rcount<10
                r1(i) =randsample(SearchAgents_no,1,true,Pr);
                rcount=rcount+1;
            else
                break;
            end
        end
        rcount=0;
  
        while first == r2(i) || r1(i) == r2(i)
            if rcount<10
                r2(i) =randi(SearchAgents_no+nA);
                rcount=rcount+1;
            else
                break;
            end
        end
        rcount=0;
        V(:, first) = X(:, first) + FFFF(1) .* (X(:, r1(i)) - PA(:, r2(i)));
        for j = 1 : dim
            if V(j, first) < lb
                V(j, first) = 0.5 * (lb + X(j, first));
            end
            if V(j, first) > ub
                V(j, first) = 0.5 * (ub + X(j, first));
            end
        end
        jrand = randi(dim);
        for j = 1 : dim
            if rand <= CCCC(1) || j == jrand
                U(j, first) = V(j, first);
            else
                U(j, first) = X(j, first);
            end
        end
    end
    fu=feval(fobj,U,fnum);
    fu=fu-fnum*100;
    counteval = counteval + SearchAgents_no;
    bR=fu<fitness;
    for i=1:H
        if sum(r==i)==0
            SR(i)=0;
        else
            SR(i)=sum(bR' & (r==i))/sum(r==i);
        end
    end
    % Selection
    nS=sum(bR);
    A= [A X(:, bR)];
    nA= nA + nS;
    S_CR	= CR(bR);
    S_F		= F(bR);
    S_df	= abs(fu(bR) - fitness(bR));
    X(:, bR)		= U(:, bR);
    fitness(bR)	= fu(bR);
    % Update MCR and MF
    if nS > 0 
        w = S_df(1 : nS) ./ sum(S_df(1 : nS));
        if all(S_CR(1 : nS) == 0)
            MCR(iM) = -1;
        elseif MCR(iM) ~= -1
            MCR(iM) = (sum(w .* S_CR(1 : nS) .* S_CR(1 : nS)) / sum(w .* S_CR(1 : nS))+MCR(iM))/2;   
        end
        MF(iM) = (sum(w .* S_F(1 : nS) .* S_F(1 : nS)) / sum(w .* S_F(1 : nS))+MF(iM))/2;  
        SR(iM)=max(SR);
        iM = mod(iM, H-1) + 1;
    end
    if SR==0
        PR=ones(1,H)/H;
    else
        PR=SR./sum(SR);
    end
    % Sort
    [fitness, fidx] = sort(fitness);
    X = X(:, fidx);
    % Update NP and population
    SearchAgents_no = round(NPinit - (NPinit - NPmin) * counteval / FEsmax);
    fitness = fitness(1 : SearchAgents_no);
    X = X(:, 1 : SearchAgents_no);
    U = U(:, 1 : SearchAgents_no);
    Asize = round(Arate*SearchAgents_no);
    if nA>Asize
        A=A(:,(nA-Asize+1):nA);
        nA=Asize;
    end
    Convergence_curve(1,countiter)=fitness(1);
    countiter = countiter + 1;
end
Global_score = fitness(1);
Global_pos = X(:, 1);
end


function x= cauchyinv(p, varargin)
% USAGE:       x= cauchyinv(p, a, b)
% Inverse of the Cauchy cumulative distribution function (cdf), x= a + b*tan(pi*(p-0.5)).
% ARGUMENTS:
% p (0<=p<=1) might be of any dimension.
% a (default value: 0.0) must be scalars or size(p).
% b (b>0, default value: 1.0) must be scalars or size(p).
% EXAMPLE:
% p= 0:0.01:1;
% plot(cauchyinv(p), p);
% SEE ALSO:    cauchycdf, cauchyfit, cauchypdf, cauchyrnd.
% Copyright (C) Peder Axensten <peder at axensten dot se>
% HISTORY:
% Version 1.0, 2006-07-10.
% Version 1.1, 2006-07-26.
% - Added cauchyfit to the cauchy package. 
% Version 1.2, 2006-07-31:
% - cauchyinv(0, ...) returned a large negative number but should be -Inf. 
% - Size comparison in argument check didn't work. 
% - Various other improvements to check list. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Default values
	a=	0.0;
	b=	1.0;
	% Check the arguments
	if(nargin >= 2)
		a=	varargin{1};
		if(nargin == 3)
			b=			varargin{2};
			b(b <= 0)=	NaN;	% Make NaN of out of range values.
		end
	end
	if((nargin < 1) || (nargin > 3))
		error('At least one argument, at most three!');
    end
	p(p < 0 | 1 < p)=	NaN;
	% Calculate
	x=			a + b.*tan(pi*(p-0.5));
	% Extreme values. 
	if(numel(p) == 1), 	p= repmat(p, size(x));		end
	x(p == 0)=	-Inf;
	x(p == 1)=	Inf;
end


function r= cauchyrnd(varargin)
% % % % %Generate random numbers from the Cauchy distribution, r= a + b*tan(pi*(rand(n)-0.5)).
% % % % Chy = cauchyrnd(0, 0.1, SearchAgents_no + 10);%(SearchAgents_no + 10)*(SearchAgents_no + 10) rand numbers
% USAGE:       r= cauchyrnd(a, b, n, ...)
% Generate random numbers from the Cauchy distribution, r= a + b*tan(pi*(rand(n)-0.5)).
% ARGUMENTS:
% a (default value: 0.0) must be scalars or size(x).
% b (b>0, default value: 1.0) must be scalars or size(x).
% n and onwards (default value: 1) specifies the dimension of the output.
% EXAMPLE:
% r= cauchyrnd(0, 1, 10); % A 10 by 10 array of random values, Cauchy distributed.
% SEE ALSO:    cauchycdf, cauchyfit, cauchyinv, cauchypdf.
% Copyright (C) Peder Axensten <peder at axensten dot se>
% HISTORY:
% Version 1.0, 2006-07-10.
% Version 1.1, 2006-07-26.
% - Added cauchyfit to the cauchy package.
% Version 1.2, 2006-07-31:
% - cauchyinv(0, ...) returned a large negative number but should be -Inf.
% - Size comparison in argument check didn't work.
% - Various other improvements to check list.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Default values
a=	0.0;
b=	1.0;
n=	1;
% Check the arguments
if(nargin >= 1)
    a=	varargin{1};
    if(nargin >= 2)
        b=			varargin{2};
        b(b <= 0)=	NaN;	% Make NaN of out of range values.
        if(nargin >= 3)
            n=	[varargin{3:end}];		
        end
    end
end
% Generate
r=	cauchyinv(rand(n), a, b);
end


